EPJ manuscript No. 

(will be inserted by the editor) 



Communities recognition in the Chesapeake Bay ecosystem by 
dynamical clustering algorithms based on different oscillators 
systems 

Alessandro Pluchino 1 , Andrea Rapisarda 1 and Vito Latora 1 

Dipartimento di Fisica e Astronomia, Universita di Catania, 

and INFN sezione di Catania, Via S. Sofia 64, 1-95123 Catania, Italy 

Received: date / Revised version: date 

Abstract. We have recently introduced [T][2] an efficient method for the detection and identification of 
modules in complex networks, based on the de-synchronization properties (dynamical clustering) of phase 
oscillators. In this paper we apply the dynamical clustering tecnique to the identification of communities of 
marine organisms living in the Chesapeake Bay food web. We show that our algorithm is able to perform 
a very reliable classification of the real communities existing in this ecosystem by using different kinds 
of dynamical oscillators. We compare also our results with those of other methods for the detection of 
community structures in complex networks. 
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1 Introduction 

Complexity theory and associated methodologies are trans- 
forming ecological research, providing new perspectives on 
old questions as well as raising many new ones. Patterns 
and processes resulting from interactions between individ- 
uals, populations, species and communities in landscapes 
are the core topic of ecology. These interactions form com- 
plex networks, which are the subject of intense research 
in complexity theory, informatics and statistical mechan- 
ics. This research has shown that complex natural net- 
works often share common structures such as loops, trees 
and clusters, which contribute to widespread processes in- 
cluding feedback, non-linear dynamics, criticality and self- 
organization. 

In ecosystems, and in particular in food webs, these struc- 
tures have strong implications for their stability and dy- 
namics. Actually, a food web constitutes a special descrip- 
tion of a biological community with focus on trophic inter- 
actions between consumers and resources [3]. Food webs 
are deeply interrelated with ecosystem processes and func- 
tioning since the trophic interactions represent the trans- 
fer rates of energy and matter within the ecosystem. In 
particular it is known that trophic webs are not randomly 
assembled, but are the result of the interaction of different 
cohesive subgroups (modules or community structures). 
Therefore, identifying the tightly connected groups within 
these networks is an important tool for understanding the 
main energy flows of the networks itself, as well as for 
defining a hierarchy of nodes and connections within a 
complex structure. 




Fig. 1. An example of a network made of four communities 
or modules, defined as subsets of nodes within which the net- 
work connections (links) are dense, but between which they are 
sparse. 



For practical purposes, modules can be defined as subsets 
of network nodes within which the connections are dense, 
but between which they are sparse (see FigJT]). In the last 
years many efficient heuristic methods have been proposed 
to investigate the presence of these structures in com- 
plex networks, and their performances have been tested 
on both real and computer generated networks with a 
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known subdivision in different communities 4.5 . We have 
recentfy presented a dynamical clustering (DC) algorithm 
for modules identification based on the de-synchronization 
properties of a given dynamical system associated to the 
network [HI3- It combines topological and dynamical in- 
formation in order to recognize modular structures with 
a precision and a computational cost (0(N 2 ) on a sparse 
graph) competitive with the majority of the other tech- 
niques. In this paper we apply the DC algorithm to a well- 
known food web of marine organisms living in the Chesa- 
peake Bay, situated on the Atlantic coast of the United 
States (see FigEJ). We implement our algorithm by using 
several dynamical systems and we compare the results of 
the simulations among them and also with those obtained 
with other methods. 

2 Dynamics of weighted networks of coupled 
oscillators 

The DC algorithm is based upon the well-known phe- 
nomenon of synchronization of coupled phase oscillators 
[6] , each one associated to a node of a given network, and 
interacting through the edges or links of the graph. In 
Ref.[7] it has been shown that an enhancement in the ca- 
pability of synchronization can be achieved by using the 
information contained in the overall topology of the net- 
work. This can be realized in practice through a weighting 
procedure wich associates a load to each link of the net- 
work. The load Uj of the link connecting nodes i and j 
can be quantified by the so called edge betweenness, i.e. 
the fraction of shortest paths that are making use of that 
link. Within this assumption, the dynamics of a network 
of N coupled oscillators {x,- is described by the 
following set of first order differential equations: 

±i = F(x t ) - - ° V 1° H[x t - xj] , (1) 

£~>j£Ki '../ jeK . 

where F = F(x) governs the dynamics of each individual 
oscillator, H = H(x) is the coupling function, a is the 
overall coupling strength and Ki is the set of neighbors 
of node i th . Notice that the loads {hj} have been raised 
to a power a, where a is a real tunable parameter which 
regulates the dynamical clustering process. 
In Ref.pQ we showed that, for a given dynamical system 
F(x;) and for a given value of the coupling strength c, 
if the system starts in a perfectly synchronized state at 
a = and a is let to slowly decrease in time from to 
— oo, the links with the higher load will be weighted less 
and less with respect to the other links, thus inducing 
a progressive desynchronization (dynamical clustering) of 
the system in a hierarchy of clusters of oscillators cor- 
responding to different configurations of modules for the 
network considered. In order to select which one of these 
configurations is the best one as a function of a(t), we 
decided to look to local or global maxima of the modular- 
ity Q. The latter simply compares the fraction of edges 
within n c arbitrary communities (intra-community links) 




Fig. 2. Geographic position of the Chesapeake Bay ecosystem, 
situated on the Atlantic coast of the United States. 



of a given network with the expected fraction of such edges 
in a random network, which does not exhibits community 
structures [8j. Actually, it is possible to define a n c x n c 
size matrix e where the elements represent the fraction 
of total links starting at a node in partition i and ending 
at a node in partition j. Then, the sum of any row (or col- 
umn) of e, ai = eij corresponds to the fraction of links 
connected to i. If the communities were allocated without 
any regard to the underlying structure, the expected num- 
ber of intra-community links would be just x a^. On the 
other hand, we know that the fraction of links exclusively 
within a partition is en. So, we can compare the two di- 
rectly and sum over all the partitions in the graph: 

Q^Y.^-°Z) ( 2 ) 

i 

Obviously it makes sense to look for high values of Q. In 
fact, if we take the whole network as a single community, 
we get Q — 0, while values approaching Q = 1 indicate 
strong community structure; on the other hand, for a ran- 
dom network we get again Q = 0. It is important to notice 
that the expression ^ is not normalized, so that Q can- 
not reach in practice the value 1. For networks with an 
appreciable subdivision in classes, Q usually falls in the 
range from about 0.2 to 0.7. 

In Refs.[l][2] we applied the DC algorithm to several real 
and trial networks, using as dynamical systems the Opin- 
ion Changing Rate and the Rossler ones and adopting 
modularity Q to choose the best subdivision for a given 
network. In the next sections, by using those and other 
dynamical systems - like the Kuramoto's one -, we will 
adopt again the modularity approach in order to explore 
the complex modular structure of the trophic relationships 
among organisms living in the Chesapeake Bay. 
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Fig. 3. Chesapeake Bay food web. The nodes represent the 
33 most important taxa, interconnected by 71 links (preda- 
tory interactions) and grouped in two main communities: the 
Benthic organisms (full circles) and the Pelagic organisms (full 
squares) . Only four taxa (open circles) are undetermined. The 
modularity of such a natural subdivision is Q = 0.337, if one 
considers nodes 17 and 28 as belonging to the Benthic commu- 
nity and nodes 6 and 11 to the Pelagic one. 

3 Dynamical clustering analysis of the 
Chesapeake Bay food web 

The Chesapeake Bay watershed is a large and complex 
ecosystem of the U.S. Atlantic Coast, made up of smaller 
subsystems including forests, streams and marshes. Ecosys- 
tems work through the plants and animals that live in 
them. In a healthy ecosystem, plants and animals can ben- 
efit each other in a cycle of energy. Plants use solar energy 
to grow, transforming nutrients from the decay and waste 
of other living things. Animals eat the plants and recy- 
cle the nutrients, through their wastes and by their death 
and decay, for the use of other living things. The same 
process occurs on the land, in terrestrial ecosystems, and 
in the water, in aquatic ecosystems. Ecosystems continue 
to thrive when the energy from the nutrients in this cycle 
is not wasted or lost, but is stored and recycled. 
The Chesapeake Bay (CB) ecosystem was originally stud- 
ied by Baird and Ulanowicz [9] , who carefully investigated 
the trophic relationships, i.e. the predatory interactions, 
between the 33 most important taxa (i.e. species or groups 
of species), which can be represented by the nodes of a 
food network with 71 links (see Fig[3]). Baird and Ulanow- 
icz compiled the matrix of these relationships specifying 
the percentage of carbon assimilated in each interaction, 
thus the network should be directed (A feeds on B but 
the opposite is not true), and valued (due to the differ- 



ent percentages of carbon exchanged in the interactions). 
However, as usually done in many papers [lOllllj . we will 
consider it non-directed and non-valued. 
In Fig|3] the two main communities of organisms classi- 
fied by Baird and Ulanowicz are visible: the Benthic ones 
(full circles), which live near the bottom of the bay, and 
the Pelagic ones (full squares), which live near the sur- 
face or at middle depths) . Only four species (open circles) 
are undetermined. We will verify in the following if the 
topological information about the predatory interactions, 
compressed in the edge betweennesses of the links of the 
CB food web (and stored in the load matrix kj , calculated 
once forever for this network), is sufficient to identify, by 
means of the dynamical clustering algorithm, these two 
main a — priori communities, or similar configurations 
preserving in some way the distinction between Benthic 
and Pelagic organisms. 

Before going on, it is important to immediately calcu- 
late, by using the definition given in the previous section, 
the modularity Q of the natural subdivision (Benthic vs 
Pelagic) shown in Fig[3l Actually, there are several pos- 
sible configurations depending on the arrangement of the 
nodes correspnding to undetermined species. If we con- 
sider (as it seems more natural looking at the network 
itself) nodes 17 and 28 as belonging to the Benthic com- 
munity and nodes 6 and 11 to the Pelagic community, the 
resulting modularity is Q = 0.337, a high value which in- 
dicates a good subdivision. On the other hand, attribut- 
ing node 11 to the Benthic community and node 28 to 
the Pelagic one would produce a lower modularity score 
Q = 0.283. Considering nodes 11 and 28 as a third sepa- 
rated community would give back a modularity Q = 0.321. 
Therefore, assuming the configuration with two commu- 
nities and Qref = 0.337 as our reference configuration, in 
the next sections we apply our DC algorithm to the Chesa- 
peake Bay food web, using different oscillators' systems, in 
order to compare the resulting best configurations among 
them and with the one chosen as reference. 



3.1 Dynamical clustering with a system of Rossler 
oscillators 

The dynamics of a system of N identical (three-dimen- 
sional) chaotic Rossler oscillators, defined over the nodes 
of a given network, is ruled by Eq. |T]) , with X; = (xt , yi, Zi) , 
F(xj) = [-uiyi - Zi,u>Xi + 0.165^,0.2 + Zi(xi - 10)] and 
H(x) = [a;, 0,0] (thus the coupling acts only on the x 
variable). In other words we have the following equations 
of motion [2]: 

±i = -uiy, -Zi- a l % ( x * " x i) 

y, = uxi + Q.WByi (3) 

^ = 0.2 + ^(^-10) i = l,...,N. 

Here a; is a common natural frequency associated at each 
oscillator that, without loss of generality, we put equal to 
1.0. As previously seen, the load matrix Zy (the matrix of 
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Fig. 4. A typical run of the DC algorithm for the Chesapeake 
Bay food web: time-evolution of both the Rossler's phases (top 
panel) and the corresponding modularity (bottom panel) as a 
function of a(i). 
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Fig. 5. A typical run of the DC algorithm for the Chesapeake 
Bay food web: time-evolution of both the Kuramoto's frequen- 
cies (top panel) and the corresponding modularity (bottom 
panel) as a function of a(t). 



Table 1. For the Rossler case, clusters configuration with the 
best modularity score Qbest = 0.43 at ot^est = —1.62, indicated 
with an arrow in the bottom panel of Fig|4] 



cluster 


nodes 


cluster 1 (10 nodes) 


3,14,15,16,18,25,26,27,28,29 


cluster 2 (3 nodes) 


4,17,19 


cluster 3 (1 nodes) 


30 


cluster 4 (3 nodes) 


22,31,32 


cluster 5 (14 nodes) 


1,2,7,8,9,10,11,12,13,20,21,23,24,33 


cluster 6 (2 nodes) 


5,6 



the edge betweennesses) is calculated once forever for the 
chosen network (in this case the CB food web). 
In order to evaluate the degree of synchronization of the 
Rossler system (|3.1j) one has to calculate the order param- 

eter W = <i| £^ e**«<*>|) tj where *,(f) = arctan[^} 
indicates the istantaneous phase of the i-th oscillator and 
{■■■)t stays for a time average. If all the oscillators rotate 
independently, no clusters exist and we have & ~ . On 
the contrary, if their motions are synchronized in phase, 
only one cluster exists and we obtain <F ~ 1. Once a net- 
work is fixed, the first task is to find the value of the cou- 
pling parameter a providing a fully synchronized starting 
state for the Rossler oscillators at a — (i.e. at t = 0). 
Then, one can let a to decrease in time and study the 
dynamical clustering process acting on the istantaneous 
phases <Pi(t)'s of the oscillators. We call "cluster" a group 
of contiguous phases in the <£'s interval (usually [—3,3]) 
separated by a distance of more than 0.02 units. For each 
value of a a different configuration of clusters (correspond- 
ing to a given network subdivision) will appear and one 
has to calculate the corresponding modularity and select 
the configuration with the best modularity score. 
In FigHwe show the result of a typical event of the Rossler 
DC algorithm for a value of the interaction strenght a = 1 



(such that the system would lie in its synchronized phase 
for a = 0). The clusters evolution (top panel) and the 
corresponding modularity Q(t) (bottom panel) are plot- 
ted as a function of a (note that, in the top panel, the 
average istantaneous phase of the system has been sub- 
tracted from the istantaneous phases of the oscillators in 
order to have a symmetric plot). The system starts in 
a fully synchronized state (xi(0) — yi(0) — Zj(0) = 
Vi) at a s tart — and evolves through decreasing val- 
ues of a(t) (with a decrement 5a = 0.0008), up to the 
value a enc i = —2. Even if the system strongly oscillates 
during the desynchronization process, clusters' configura- 
tions (i.e. community structures of the underlying net- 
work) with large values of modularity appear. The de- 
tailed configuration with the highest modularity peak (see 
the arrow in the bottom panel) is reported in Table [T] 
It consists of 6 clusters with a Qbest = 0.43, obtained 
for ctbest — —1.62, and it is quite consistent with the 
distinction between pelagic organisms (clusters n.1,2 and 
4) and benthic organisms (clusters n.3,5 and 6). Further- 
more, if compared with the reference configuration, where 
Qref = 0.337, the configuration we found here seems ev- 
idently (having a higher modularity) to better reflect the 
underlying structure emerging from the global information 
stored in the edge betweennesses of the food web. 



3.2 Dynamical clustering with the Kuramoto model 

The Kuramoto model describes a population of N peri- 
odic oscillators having natural frequencies uii and coupled 
through the sine of their phase differences [12] . It is simple 
enough to be analytically solvable, still retaining the basic 
principles to produce a rich variety of dynamical regimes 
and synchronization patterns |13pi4j . 
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Table 2. Fot the Kuramoto case, clusters configuration with 
the best modularity score Qbest = 0.404 at ctb es t = —3.50, 
indicated with an arrow in the bottom panel of Fig(5] 



cluster 


nodes 


cluster 1 (2 nodes) 


4,17 


cluster 2 (2 nodes) 


9,10 


cluster 3 (11 nodes) 


3,14,15,16,18,19,25,26,27,28,29 


cluster 4 (14 nodes) 


1,2,7,8,11,12,13,20,21,22,23,24,30,33 


cluster 5 (2 nodes) 


5,6 


cluster 6 (2 nodes) 


31,32 



The dynamics of the model is given by 
K N 

6i(t) =Ui + jfJ2 sin ^" ~ *<) * = " • • ' N W 
3=1 

where 6i(t) is the phase (angle) of the ith oscillator at 
time t, while is its intrinsic frequency randomly drawn 
from a symmetric, unimodal distribution g(oS) with a first 
moment u>o (typically a Gaussian distribution or a uniform 
one) . These natural frequencies u>i are constant and time- 
independent. The sum in the above equation is running 
over all the oscillators so that this is an example of a 
globally coupled system. The most interesting feature of 
the model is that, despite the difference in the natural 
frequencies of the oscillators, it exhibits a spontaneous 
transition from incoherence to collective synchronization 
beyond a certain threshold K c of the coupling strength 
K [14] . For small (positive) values of K, each oscillator 
tends to run independently with its own frequency, while 
for large values of K, the coupling tends to synchronize 
(in phase and frequency) the oscillator with all the others. 
In the Kuramoto model, it is possible to define a complex 
mean field order parameter as 

3=1 

where the magnitude < r(t) < 1 is a measure of the co- 
herence of the population and \P(t) is the average phase. 
In other words, as N approaches infinity, the magnitude 
roothe average istantaneous phase of the system has been 
subtracted from the istantaneous phases of the oscillators 
in order to have a symmetric plot) of the complex mean 
field after a transient time should be zero in the incoherent 
state with K < K c and different from zero in the coher- 
ent state with K > K c . Actually, as K increases beyond 
K c , more and more oscillators will be recruited toward 
the mean phase ^(i) and is expected to continuously 
increase from zero to one. 

In order to use Kuramoto model as dynamical system 
in the context of our dynamical clustering algorithm, we 
immediately note that, at variance with the Rossler sys- 
tem, Eqs|4| are already very similar to Eqs[Tl providing 
that one considers jq = 6*;, F(x^) — uji, a — K, Uj = 1 
and H(x) = sm(0j — In other words Kuramoto equa- 
tions already have a coupling term, containing a nonlinear 



function that becomes approximatively linear in the syn- 
chronization manifold (where 8j ~ This makes Ku- 
ramoto model particularly suitable for our purpose. In fact 
in this case, after having chosen a given network and hav- 
ing associated an oscillator to each node, we can directly 
apply the weighting procedure of section 2 (see EqQ} to 
Eqs[4]without adding further coupling terms, thus simply 
obtaining: 

6i(t) =Wi+ ^ l % sin (^ ~ B *) i = 1, • • • , JV 

(6) 

where as usual, Uj is the load of the link connecting nodes 
i and j in the chosen network, Ki the set of neighbors of 
node i th and a is a real tunable parameter. 
As in the previous section, our task is now to test the sen- 
sitivity of the Kuramoto dynamical clustering algorithm 
on the CB food network. First of all, we have to fix the 
coupling parameter K of Eq[6] in order to obtain a fully 
synchronized state of the system for a ~ 0. We found that 
for K > 5 such a state is guaranteed, thus we will reason- 
ably set K = 10. In our simulations of the Kuramoto 
system we will always use as initial conditions uniform 
distributions for both the OiS (in the interval [— it, ir]) and 
a>i's (in the interval [—2, 2]). We remind that the latter are 
constant in time. At variance with the Rossler case, we are 
now interested to the istantaneous frequencies 9i(t) of the 
oscillators (which at t = coincide with the natural fre- 
quencies uii). Again, the average istantaneous frequency of 
the system will be subtracted from the istantaneous fre- 
quencies of the oscillators in order to have a symmetric 
plot. 

In Fig[5] we show a typical run of the Kuramoto DC al- 
gorithm for the CB food network. As before, the simu- 
lation starts from a star t = then a decreases in time 
with a given step (5a = 0.0008): one sees that in a few 
steps the system suddenly synchronizes (due to the high 
value of K) then slowly relaxes producing a progressive 
desynchronization that is, again, very oscillating in time; 
for each value of a(t) the istantaneous clusters' configu- 
ration of frequencies is identified (being the definition of 
'cluster' the same than in the previuos section), and the 
correspondent modularity calculated. The detailed config- 
uration with the highest modularity peak (see the arrow 
in the bottom panel) is reported in Tabled It consists of 6 
clusters with a Qbest — 0.404, obtained for a& e st = —3.50: 
even if this value of modularity is less than in the Rossler 
case, on the other hand it is greater than Q re f = 0.337 and 
in any case is again quite consistent with the distinction 
between pelagic organisms and the benthic ones. 

3.3 Dynamical clustering with the Opinion Changing 
Rate model 

As a last example of application of the dynamical clus- 
tering algorithm to the Chesapeake Bay food web, let us 
to consider as dynamical system the so called Opinion 
Changing Rate (OCR) model PQ. It was originally in- 
troduced in Ref.[15] as a modification of the Kuramoto 
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model, in order to study how the personal inclination to 
change, randomly distributed in a group of individuals, 
can affect the opinion dynamics of the group itself. The 
dynamics of a system of N fully coupled individuals (os- 
cillators) is governed by the following set of differential 
equations: 



Xi(t) =Ui + -^Yl P sin ( x r 



. Xi ) e -P\*i-*i\ i = l,.,.,N 



(7) 

where Xi(t) is the opinion of the ith agent at time t. 
Here the opinions have a very general meaning and can 
be usefully represented by means of unlimited real num- 
bers Xi s] — oo + oo[ Vi = 1, ...,N. Opinions interact 
by means of the coupling term, where a is the coupling 
strength and the exponential factor, tuned by the param- 
eter /?, ensures that opinions will not influence each other 
any longer when the reciprocal distance exceeds a certain 
threshold. This is perhaps the most remarkable feature 
of the OCR model, since it allows the system to reach 
an asymptotic stationary state where the configuration of 
opinions does not vary any longer. The parameter /3 ap- 
pears also as a factor of the sine in the coupling term and 
simply rescales the range of the coupling strenght. We typ- 
ically adopted the value /3=3, which ensures a consistent 
behavior of the exponential decay. Finally, the u>i's - cor- 
responding to the natural frequencies of the oscillators in 
the Kuramoto model - represent here the so-called natural 
opinion changing rates, i.e. the intrinsic inclinations of the 
agents to change their opinions. For this reason we called 
this model the Opinion Changing Rate (OCR) model |15j . 
The values lu^s, which do not depend on time, are uni- 
formly distributed in the range [—0.5, 0.5] and represent 
also the initial conditions for the istantaneous frequencies 

Xi(t)'s. 

Many numerical simulations were performed in Ref.[15j 
starting from a uniform distributions of the initial opin- 
ions Xi(t = 0) in the range [—1,1]. As a function of the 
coupling strength a, a transition was observed from an 
incoherent phase (for a < a c ), in which people tend to 
preserve different opinions and different frequencies ac- 
cording to their natural changing rates u>i, to a partially 
synchronized phase, where people share a small number 
of opinions, and, finally, to a fully synchronized one (for 
a » a c ) in which all the people share the same opin- 
ion changing with the same rate. In order to measure 
the degree of synchronization of the system, it can be 
adopted an order parameter related to the standard devia- 
tion of the istantaneous frequencies and defined as R(t) = 

1 — J jj ^2f = i(ij(t) — X(t)) 2 , where X(t) is the average 

over all individuals of x'j(t). From such a definition it fol- 
lows that R = 1 in the fully synchronized phase and R < 1 
in the incoherent or partially synchronized phase. 



3.3.1 Standard Opinion Changing Rate model 

In order to utilize the OCR model as dynamical system for 
recovering community structures in the CB food network 
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Fig. 6. A typical run of the DC algorithm for the Chesapeake 
Bay food web: time-evolution of both the OCR's frequencies 
(top panel) and the corresponding modularity (bottom panel) 
as a function of a(t). 

Table 3. For the OCR model, clusters configuration with the 
best modularity score Qbest = 0.404 at ctbest ~ —5.89, indi- 
cated with an arrow in the bottom panel of Fig(6] 



cluster 



nodes 



cluster 1 (14 nodes) 
cluster 2 (2 nodes) 
cluster 3 (11 nodes) 
cluster 4 (2 nodes) 
cluster 5 (2 nodes) 
cluster 6 (2 nodes) 



1,2,7,8,11,12,13,20,21,22,23,24,30,33 
31,32 

3,14,15,16,18,19,25,26,27,28,29 

5,6 

9,10 

4,17 



we put togheter EqfT]and EqJTl thus obtaining 



x\(t) 



E 



jZKi ij j eK . 



where i = 1, . . . , N, a is the usual real tunable parameter 
and Ki is the set of neighbors of node i . As in the case 
of Kuramoto model, we do not need any further coupling 
term in Eqs. [SJ since such a term is already present in 
the OCR model. We will follow now the time evolution of 
the istantaneous frequencies x\(t) starting, as usual, from 
a completely synchronized state at a ~ 0: more precisely, 
as in the Kuramoto case, the initial frequency distribution 
is uniform inside the interval [—0.5,0.5] (since it coincides 
with the LUi's distribution) but we chose a — 5.0, a cou- 
pling that ensures a rapid synchronization for t > 0. Then 
we let a to decrease in time with a step 5a = 0.01. We are 
confident that, at variance with the previous examples, the 
exponential factor in the coupling term could trigger the 
aggregation of frequencies in stable clusters corresponding 
to community configurations of the network. 
In Fig[6]we plot the OCR frequencies' time evolution for 
the Chesapeake Bay food web, together with the respec- 
tive modularity. One can immediately appreciate the sta- 
bility of the desynchronization process, which produces 
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(due to the exponential factor in the coupling term) a 
branching sequence of metastable plateaux correspond- 
ing to different clusters configurations. The best one is 
obtained for —5.9 > a^est > —8.3 with a modularity 
Qbest = 0.404 and corresponds exactly to the same best 
configuration (six clusters) of the Kuramoto case, see Ta- 
blej3] But the stability of the synchronized manifold makes 
useless, in this case, to adopt lower values of 6a, thus mak- 
ing also the simulation less expensive in terms of compu- 
tational cost. 



3.3.2 Opinion Changing Rate model with HK dynamics 

In order to further improve the performances of the OCR 
system, in Ref. pQ we tought to modify the natural opinion 
changing rates w's following a suggestion from the Hegsel- 
mann and Krause (HK) model of opinion formation. The 
HK model [H] is based on the concept of bounded con- 
fidence, i.e. on the presence of a parameter e, called con- 
fidence bound, which expresses the compatibility among 
the agents in the opinion space. If the opinions of two 
agents i and j differ by less than e, their positions are close 
enough to allow for a discussion, which eventually leads 
to a change in their opinions, otherwise the two agents do 
not interact with each other. In our case, the opinion space 
is one-dimensional, being usually Xi(t) G [—1,1]. Thus we 
fixed a small value for e and we let the w's to change in 
time (while in the standard OCR version the w's are time 
independent) starting from a random uniform distribution 
in the interval [—0.5, 0.5] and according to the following 
HK dynamics: at each time step a given agent, with an 
opinion Xi and a natural frequency u>i, checks how many 
of its neighbors (according to the network topology) are 
compatible, i.e. lie inside the confidence range [xi — e, Xi+e] 
in the opinion space. Next, the u>i(i) of the agent takes the 
average value of the w's of its compatible neighbors at time 
t — 1. We will refer to this new system as OCR — HK. In 
the OCR-HK system the secondary process acting on the 
w's is superimposed to the main dynamical evolution of 
Eq[5] and contributes to further stabilize the istantaneous 
frequencies Xi{t) of the agents (oscillators). 
In FigtZl the new simulation with the OCR-HK system 
are shown for the Chesapeake Bay food web. This simu- 
lation was performed again with a = 5.0 and a 5a = 0.01, 
with a confidence bound e = 0.0005. Again a branching 
desynchronization process with well defined metastable 
plateaux occurs. But in this expected, the HK 

dynamics produces a highest value of modularity, Qbest = 
0.42 (very near to that obtained with the Rossler algo- 
rithm), which is reached for —6.77 > a^est > —10.67, 
yielding a subdivision of the food web into five communi- 
ties, whose detailed structure is shown in TableU It clearly 
appears that, besides the greater value of Q, the distinc- 
tion between pelagic and benthic organisms is improved 
with respect to the standard OCR case too. In any case, 
both the standard OCR and the OCR-HK algorithm give 
a best configuration whose modularity is higher than the 
reference one Q re f = 0.337. 
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Fig. 7. A typical run of the DC algorithm for the Chesapeake 
Bay food web: time-evolution of both the OCR-HK's frequen- 
cies (top panel) and the corresponding modularity (bottom 
panel) as a function of a(t). 

Table 4. For the OCR-HK case, clusters configuration 
with the best modularity score Qbest = 0.42 at atest ~ 
—6. 77, indicated with an arrow in the bottom panel of Fig[7] 



cluster 



nodes 



cluster 1 (2 nodes) 
cluster 2 (11 nodes) 
cluster 3 (2 nodes) 
cluster 4 (14 nodes) 
cluster 5 (2 nodes) 



31,32 

3,14,15,16,18,19,25,26,27,28,29 
4,17 

1,2,7,8,9,10,11,12,13,20,21,22,23,24,30,33 
5,6 



4 Discussion and Conclusions 

Summarizing, the dynamical clustering (DC) algorithm, 
which exploits the synchronization properties of some dy- 
namical system of oscillators associated with the nodes of 
a given complex network, seems to work very well in iden- 
tifying the underlying community structure of the Chesa- 
peake Bat food web of trophic relationships. In fact we 
shown that, whatever the dynamical system we use (Rossler, 
Kuramoto, OCR or OCR-HK), the DC algorithm discov- 
ers community configurations with high values of modular- 
ity, in all the cases higher than the reference configuration 
of Fig|3] related to the main subdivision of the Chesapeake 
Bay network in Benthic and Pelagic organisms. This would 
imply also that, if we consider modularity optimization 
as a valid method to retrieve the best community struc- 
ture compatible with the information stored in the topol- 
ogy of a given network, then we should conclude that the 
rigid subdivision in Benthic and Pelagic organisms is not 
the optimal one for the Chesapeake Bay food web. Ac- 
tually, modularity method has been recently questioned 
by [17] , which showed that it is a — priori impossible to 
tell whether a module, detected through modularity opti- 
mization, is indeed a single module or a cluster of smaller 
modules, but the problem is still open. In any case it is 
worthwhile to notice that the DC algorithm produces bet- 
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ter results also if compared with other methods for detect- 
ing community structures using modularity optimization, 
which in the past were applied to the Chesapeake Bay 
(CB) food web. For example in Ref.[8], where Girvan and 
Newman presented their original iterative method, based 
in finding and removing progressively the edges with the 
largest betweenness until the network breaks up into its 
components, the best modularity value for the CB food 
web was Qgn = 0.380; and in Ref.pT], where another 
method based on the Information Centrality, was pro- 
posed and applyied to the CB food web, the best score 
for the modularity was Qic — 0.376. 
In conclusion, the DC algorithm, in particular in the ver- 
sion using the OCR-HK system (which produces large 
modularity configurations with very stable synchroniza- 
tion patterns), seems to be a very efficient method for the 
study of community structures in ecosystems and food 
webs, even if in the approximation of non-directed net- 
works. In this direction, a further improvement in the 
DC algorithm performance probably could come from the 
use of a recent generalization of the modularity approach 
which incorporate also the information contained in edge 
directions [18] , 
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